Thermoreversible Associating Polymer Networks: I. Interplay of Thermodynamics, 

Chemical Kinetics, and Polymer Physics 
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Hybrid molecular dynamics/Monte Carlo simulations are used to study melts of unentangled, ther- 
moreversibly associating supramolecular polymers. In this first of a series of papers, we describe 
and validate a model that is effective in separating the effects of thermodynamics and chemical 
kinetics on the dynamics and mechanics of these systems, and is extensible to arbitrarily nonequi- 
librium situations and nonlinear mechanical properties. We examine the model's quiescent (and 
heterogeneous) dynamics, nonequilibrium chemical dynamics, and mechanical properties. Many of 
our results may be understood in terms of the crossover from diffusion-limited to kinetically-limited 
sticky bond recombination, which both influences and is influenced by polymer physics, i. e. the 
connectivity of the parent chains. 

PACS numbers: 83.80.Kn,83.10.Rs,82.35.-x,81.05.Lg 



I. INTRODUCTION 

Flexible synthetic polymers have long been of funda- 
mental scientific interest because many of their properties 
arise from a few universal features like the topological 
connectivity, random-walk like structure, and excluded 
volume of the chain molecules. Less universal are the var- 
ious attractive, "associative" intermolecular interactions 
[l|, 0] ranging from weak dispersion forces to strong cova- 
lent chemical bonds (in chemically crosslinked systems). 
Examples include hydrogen bonding, electrostatic attrac- 
tions, and effective attractions driven by incompatibility 
with a solvent. These interactions lead to formation of 
supramolecular structures ranging from micelles to net- 
work gels. 

Associating polymers (APs) differ from simple ho- 
mopolymers in that chains contain a (typically fairly low) 
fraction of "sticky" monomers, which are different from 
the majority-species monomers. The sticky monomers 
form "sticky bonds" with each other via associative in- 
teractions weaker than permanent covalent bonds. These 
lead to formation of supramolecular aggregrates. Unlike 
the closely related "living" or "equilibrium" polymers, 
the degree of polymerization of AP "parent" chains is 
fixed (in time) by permanent covalent backbone bonds. 

The lifetime of the sticky bonds is finite. Depending 
on the nature of the associative interactions and am- 
bient conditions (e. g. temperature, concentration), the 
supramolecular topology may be practically permanent, 
in which case the system forms a "chemical gel" (i. e. a 
crosslinked rubber), or so short-lived that the system is 
indistinguishable from a simple polymer solution, melt 
or glass. Between these limits, when the topology of the 
associated supramolecular aggregates changes on a time 
scale comparable to the experiment, these systems form 
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complex fluids with fascinating dynamical and mechani- 
cal properties [|[. 

At fixed ambient conditions, the time scales for topo- 
logical changes in associating polymer systems are in 
principle set by three independent factors: the (a) ther- 
modynamics (i. e. energetic strength relative to ksT), 
(b) the "chemical kinetics" of the sticky bonds, and (c) 
the underlying non- associative polymer physics. Ther- 
modynamics set "static" quantities such as the size of 
the supramolecular aggregates and hence the position of 
the system relative to the percolative gelation transition. 
Kinetics set relaxation times through their effect on the 
rates of formation and breaking of sticky bonds. Polymer 
physics alters the dynamics through such effects as the 
random- walk- like structure and uncrossability of chains, 
which give rise to the systems' underlying Rouse or rep- 
tation dynamics The interplay of (a)-(c) allows for 
the design of materials with exquisitely tunable rheologi- 
cal response. For this reason, APs have been the focus of 
intense experimental and theoretical study over the past 
two decades; see Refs. 0, H, H, H, 0] for reviews. 

Changes in ambient conditions lead to "thermore- 
versible" property changes unique to AP systems, e. g. 
extremely sharply decreasing viscosity upon decreasing 
concentration or increasing temperature. These changes 
can be tuned (engineered), so APs have great potential 
as "smart" materials @, H, S H m which the change of 
lifetime or concentration of the sticky bonds with ambi- 
ent conditions leads to useful products. Applications in- 
clude temperature-sensitive adhesives, coatings for heat- 
sensitive materials, and generally enhanced melt process- 
ability relative to conventional polymers Q . 

Static thermodynamic properties of AP systems, in 
particular the percolative gelation transition and local 
structural changes arising from associative interactions, 
have been extensively studied. Analytic theories pro- 
vide a good understanding of homogenous systems, and 
emerging numerical techniques such as self consistent 
field theoretic simulations If Oil and reaction-ensemble 
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DPD [TO] show promise for investigating inhomogeneous 
systems. 

However, the dynamical, mechanical, and nonequilib- 
rium properties of associating polymer gels and networks 
remain poorly understood. Time dependent properties 
obviously depend on kinetics, and in addition to being 
of fundamental scientific interest, a better understanding 
of them may prove important in developing new applica- 
tions of associating polymer systems such as self healing 
materials Q. The situation becomes particularly com- 
plex when the lifetime of the sticky bonds is not long 
or short compared to the "polymeric" relaxation times; 
this regime has been studied rather extensively for linear 
equilibrium polymers (see e. g. Refs. [3 El QUE 13), 
but much less so for networks. 

Analytic and quasi-analytic approaches to AP dynam- 
ics and mechanics, e. g. 1171 R efs. 0, El El in nn m, 

HI, HI IM [H, HI SlllfTIol, have made many useful, 
experimentally verifiable predictions, including nonlinear 
behaviors such as shear thickening and strain hardening 
[13, Ull- In the general case, however, the complex inter- 
play of sticky bond thermodynamics and kinetics with 
the underlying polymer physics in these systems is al- 
most certainly beyond the reach of analytic theory. For 
the sake of tractability, theories have generally neglected 
one or more features of AP systems that are likely essen- 
tial to capturing their behavior under certain ambient 
conditions. For example, as temperature drops towards 
the glass transition, attractive, non-associative interac- 
tions, such as van der Waals forces between non-sticky 
monomers, become increasingly important [32l |. More- 
over, virtually all analytic treatments have thus far been 
restricted [33j to homogeneous AP systems; the majority 
focus on the "telechelic" case of APs with only 2 sticky 
monomers per parent chain (one on each end). We be- 
lieve that inhomogeneous AP systems are the potentially 
the most interesting and useful, e. g. because inhomo- 
geneities serve to localize sticky monomer concentration 
and network connectivity, which in turn can broaden the 
relaxation spectrum [34], [35[ . 

The above set of potentially essential features of APs is 
not treated microscopically by existing theories, but can 
be readily captured by particle-based simulations. This 
is the first of a series of simulation studies, the goal of 
which is to elucidate the separate effects of sticky bond 
thermodynamics, kinetics, and other underlying polymer 
physics on the dynamical, mechanical, and nonequilib- 
rium properties of associating polymers. 

Previous particle based simulations of AP networks fo- 
cusing on d^amicalproperties, e. g. [ItJ Refs. dEHHH 
HI, Sfl EH, S3, S, H SI, EE S3 have produced a wealth 
of interesting ideas and results. However, for various rea- 
sons outlined in Section III B[ the methods employed in 
these previous studies are not suitable for the full range of 
problems we wish to consider. A more versatile model for 
APs requires the combination of: i) realistic dynamics, 
ii) applicability to far-from-equilibrium conditions, iii) 
controllable sticky binding topology, iv) variable chem- 



ical kinetics and v) the ability to treat inhomogeneous 
systems. 

In this paper we develop and validate a model, based 
on those employed in Refs. [3, SI, El], that shows all 
these characteristics. Much of the physics of amorphous 
polymer melts and gels is independent of chemical de- 
tail [1 Sol , so properties i) and ii) are captured by us- 
ing a bead-spring model [48| for the underlying (non- 
associative) polymer physics. We capture properties (iii)— 
(v) by employing a hybrid molecular dynamics / Monte 
Carlo (MD/MC) approach with controlled (specifically, 
binary) bonding and variable chemical kinetics [l4j . The 
combination of properties (i-v) allows the model to be 
used to obtain many results unattainable with previous 
methods. A key advance is that our algorithm is par- 
allclizablc. In validating and investigating our model, 
we separate the effects of thermodynamics, kinetics, and 
polymer physics on AP network dynamics and mechanics. 
We show that sticky bonding is mappable to a "mean- 
field" two-state Arrhenius rate model, but that the ki- 
netic rate constants for SB association and dissociation 
are affected by the fact that the SMs are embedded in 
polymers. 

One of the key aspects of associating polymer net- 
works, which has rarely been examined (for networks) 
in the theoretical or simulation literature, is that sticky 
bond relaxation can be either kinctically limited (i. e. 
limited by the intrinsic rate of bonding and/or debond- 
ing) or diffusion limited (i. e. when the kinetics are so 
fast that sticky bond scission and recombination events 
become correlated because newly broken SM pairs tend 
to recombine before exploring the network cage). Most 
experimental systems seem to be kinetically limited, and 
this is the case treated by almost all published analytic 
theories for reversibly associa ting networks, including 
transient network models I22I 251. Theories for kinet- 



ically limited AP networks iHaEE H, HE H3 assume 
that the sticky bond lifetime r s b is so long compared 
to all "polymeric" relaxation times that it controls all 
important time scales for network relaxation, and there- 
fore that kinetics only affect key network relaxation times 
through their effect on r s f, (or, alternately, as suggested 
by recent experiments [Hoi . l5lj . the inverse dissociation 
rate constant k^ 1 , defined below). We show that the va- 
lidity of this assumpton is questionable, and that there 
is a wide parameter space, plausibly accessible in exper- 
iments, where it is invalid. Note that published theories 
[H [H HI HE [HE for diffusion-limited reactions in 
polymeric systems, such as those of O'Shaughnessy et. 
at, either are not immdeiately applicable or have not yet 
been applied to reversible networks. 

Our simulations explore the crossover between the ki- 
netically limited and the diffusion limited cases. We con- 
firm a key prediction of Rubinstein and Semenov 24j on 
the role of bond recombination on AP dynamics, specifi- 
cally that recombination effectively renormalizes the SB 
lifetime in systems where the sticky bonds are sufficiently 
strong. However, we show that recombination couples in- 



3 



terestingly to the diffusive-kinetic crossover in a way not 
previousiy predicted. 

We extensively examine the dynamics in quiescent, 
equilibrium systems. One of the more interesting re- 
sults is that slowing chemical kinetics increases dynam- 
ical heterogeneity [401 ] in a manner similar to increasing 
the thermodynamic strength of the sticky bonds. We 
also examine nonequilibrium 'chemical dynamics' (i. e. 
in systems where the initial sticky bond population is 
not equilibrated), and nonlinear mechanical properties. 
The studies of nonequilibrium systems and mechanical 
properties presented here are limited in number because 
this paper is intended primarily to illustrate the broad 
utility of the model. 

The organization of the rest of this paper is as follows. 
In Section [TT] we further motivate, describe and exten- 
sively validate our model. We also discuss how it differs 
from those used in previous simulations of APs. In Sec- 
tion lllll we present various results for the equilibrium and 
nonequilibrium physics of thermoreversible AP networks, 
and compare to theoretical predictions. Finally, Section 
IIVI presents a discussion and conclusions. Two Appen- 
dices include technical details of the analyses. 



II. MODELS AND METHODS 

Of particular interest for the current study, and mo- 
tivating our modelling approach, are recent experiments 
performed by the group of Stephen L. Craig [3, HO, [U, 
l57l . l58l [5^ | . These have attempted to independently vary 
thermodynamics and chemical kinetics by making sys- 
tematic changes in sticky monomer chemistry (based on 
metal-ligand interactions). Systems with similar static 
properties show dramatic differences in time-dependent 
properties that are directly associated with the different 
kinetics. This effect is quite challenging to capture ex- 
perimentally because thermodynamics and kinetics are 
highly correlated for most AP systems (i. e. stronger 
binding <-> slower kinetics Q), but relatively easy to im- 
pose in simulations. 

The advantages of including variable kinetics in the 
model are rather obvious given the above discussion. 
The advantage of imposing binary bonding (i. e. a SM 
is at all times bonded to either or 1 other SM) is 
also relevance to current experiments. Sticky monomers 
with binary bonding are considered particularly valu- 
able for making thermoplastic elastomers with enhanced 
melt processability and controllable network architecture 
[60t IrJlj . Real examples include multiple hydrogen bond- 
ing monomers such as ureidopyrimidinone (UPy), which 
have highly directional associative interactions, and are 
of a strength such that the sticky bonds they form con- 
stitute a "reversible alternative for the covalent bond" 
d, i, H [H In the model used here [H, sticky 

bonds differ from covalent bonds only by their reversibil- 
ity and strength. 



A. Hybrid MD/MC Simulation Protocol 

Our model is built on the framework of the Kremer- 
Grest bead-spring model [H] , which has been extensively 
validated and is known to cap ture the key physics of 
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[63j as well as perma- 
Each associating poly- 



mer chain is linear and contains N beads (monomers). 
Systems consist of N c h chains, so the total number of 
monomers is N c hN. Periodic boundary conditions are 
imposed in all three directions, with periods Li along di- 
rections i = x, y, and z. Values of N c h range from 700 to 
5600; the lowest values are chosen to satisfy Li > 2R ee 
(where R ee is the mean chain end-end distance) , prevent- 
ing self interactions. Remaining leading order finite size 
effects in networks scale as (A^A^) -1 / 3 [65| and should 
be ~ 2% for the systems considered here. 

All monomers have mass m and interact via the trun- 
cated and shifted Lennard- Jones (LJ) potential Ulj{v) = 
4u [(a/r) 12 - (a/r) 6 - (a/r c ) 12 + (a/r c ) 6 ], where r c is 
the cutoff radius and Ui,j{r) — for r > r c . Cova- 
lent bonds between adjacent monomers on a chain are 
modeled using the finitely extensible nonlinear elastic po- 
tential U FE ne(t) = -(l/2)(kR§)ln(l - (r/i? ) 2 ), with 
the canonical [481 ] parameter choices Ro — 1.5a and 
k = SOuq/o 2 . In this study, following the majority of 
bead-spring studies on permanently crosslinked systems 
(e. g. Refs. [H], Hll), we employ flexible chains with no 
angular potential. We express all quantities in units of 
the LJ bead diameter a, intermonomer energy uq, and 
the LJ time tlj = \Jma 2 Juq. 

All systems have monomer density p — 0.85/a 3 . We 
employ two temperatures in this study: ksT = l.Oito 
and ksT = 0.6uo- These ambient conditions both corre- 
spond to dense polymer melts far above the glass tran- 
sition temperature T g ; T g ~ 0.35un /fc/? for r c — 1.5a 
and decreases with decreasing r c [63, . This far above 
T g , melt physics is known to be dominated by the repul- 
sive part of the intermonomer interactions [J]; for con- 
venience, and following convention [48j |. we use purely 
repulsive LJ interactions with However, in- 

cluding attractive interactions by increasing r c is trivial, 
is important to realistically capture T dependent prop- 
erties, and will be done in upcoming studies. 

All simulations are performed using an enhanced ver- 
sion of the LAMMPS [79[ MD code. Newton's equations 
of motion are integrated with MD using the velocity Ver- 
let method [69[ and typical timestep St — .OItlj [z3]. A 
Langevin thermostat [7l| is used to maintain the tem- 



perature. The damping time tl 
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IOOtxj is 

larger than the value typically used (7X ons ~ Tlj) in 
bead spring studies; this reduces undesirable thermostat- 
driven effects such as alteration of stress relaxatio n by 
suppression of hydrodynamic momentum transfer [72|. 
In this study we employ two "chain lengths" . Most stud- 
ies are performed at A^ = 50, which is at or below best es- 
timates of the entanglement length 50 5, N„ < 85 [73l.f74T]. 
so the melts can be fully equilibrated by allowing chains 
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to diffuse severai R ee [75j]. N = 50 is aiso a convenient 
choice of chain length because it has been considered in 
many previous studies. To elucidate the effects of under- 
lying polymeric structure on AP physics, we also consider 
monomeric melts (N — 1), that reversibly sticky-bond 
into dimers. 

After the melts are equilibrated, we choose a fraction 
c s t of the monomers to be "sticky". For the N = 1 
systems these are chosen randomly. For N = 50 systems, 
SMs are placed uniformly along chains: at both chain 
ends and also at internal monomers iN/ (Nc st — 1), where 
i = 1, 2, Nc st — 2. In this study we use c st = 0.08 
(which is comparable to typical experimental values, e. 
g. [H, |Zf|, so for N = 50 the SMs are the 1st, 17th, 
34th, and 50th monomers in each chain. However, any 
SM placement can be used, and studies of the effects 
of altering SM placement at fixed c s t are underway; the 
effects of chemical disorder are known to be significant 
for stress relaxation [29|, HH . 

Sticky monomers are identical to regular monomers, 
except that they form reversible "sticky" bonds. Fig- 
ure [T] illustrates the potential energy between sticky 
monomers as a function of their separation r. SMs (like 
all monomers) always interact via Lennard Jones interac- 
tions, whether bonded or not. Bonded SMs additionally 
interact via the potential U s b(r, h): 

U sb {r, h) = Ufene ( r ) - Ufene(i"o) ~ h, (1) 

which is based on the standard covalent FENE poten- 
tial. Here tq represents the equilibrium FENE bond 
length; ro ~ .96a, i. e. the minimum of the poten- 
tial Ulj{t) + U fene{t)- The only difference between 
the sticky and covalent bond potentials is thus an r- 
independent, tunable offset. The same bonding potential 
was used in Huang et. a/.'s studies of equilibrium poly- 
mers [13, Iz3 and a very similar potential was used in 
Baljon et. al. 's studies of telechelic associating networks 
[46l [78[ . However, our method has several important dif- 
ferences from those of Refs. [1J, |46j, [77J fsee Section lHBD . 
so we explain it in detail below. 

The potential U s b(r) has several other important fea- 
tures, h represents the sticky binding energy; for h = a 
sticky bond can be formed between two monomers sep- 
arated by ro with no change in energy. The associated 
force F s b(r) = —dU s b/dr, however, is independent of h. 
Adjusting h is thus a nearly pure way of adjusting the 
thermodynamics of the sticky bonds without directly al- 
tering their "chemical kinetics" , i. e. the rates of forma- 
tion and dissociation of sticky bonds. 

Formation and breaking of sticky bonds is performed 
using Metropolis Monte Carlo [69(. The change in en- 
ergy required to form a sticky bond between an unbonded 
pair of SMs is just AE(r, h) = U s b(r, h), and the energy 
change to break a sticky bond is AE(r, h) = —U s b(r, h). 
These are the only MC "moves" used, and the accep- 
tance probability of the moves is set by AE(r, h)/kBT. 
The MC moves are strictly "topological" . All spatial mo- 
tion of bonded SMs is governed by the sticky bond force 
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FIG. 1: Sticky monomer interaction potential. Same as in 
Refs. 0, EE]- Differences in SB formation/breaking rules are 
noted in the text. 



F s b(r), along with the other forces from Lennard Jones 
and covalent FENE interactions, which are all integrated 
using MD. One potential difficulty is that only bonded 
SMs "feel" the force F s b(r), so formation/breaking of 
SMs creates temporal force discontinuities. However, as 
will be shown below, this does not seem to cause any 
spurious behavior. 

The system sizes and time scales studied here require 
simulation on parallel computers; 8 to 64 processors are 
used in a typical simulation. While MD parallelizes very 
well it has long been noted that MC [8(| is very 
difficult to parallelize. We therefore perform hybrid par- 
allel MD / serial MC simulations. MC moves are per- 
formed once every ro in Lennard Jones time units; ro 
is the MC "timestep". The MD simulation is paused 
while the parallel-distributed lists of sticky bonds and 
SM coordinates are gathered onto one processor. For ef- 
ficiency, Verlet-style pair neighbor lists of SMs are used 
and of open SM pairs, only those within r < Rq are 
considered for SB formation. After the SB list is up- 
dated, it is distributed back to all processors and the 
MD simulation resumes. Great care was taken to opti- 
mize the MC algorithm to minimize "dead" time on the 
other processors, but reasonable parallel performance re- 
quires ro 3> St. In this paper, except where otherwise 
noted, we use r = t^j; see also Section HlDI 

As mentioned above, current experimental trends fa- 
vor binary-bonding sticky monomers, We impose binary 
bonding through a simple restriction on the Monte Carlo 
routine; sticky bond formation is attempted only for pairs 
of unbonded SMs. The 1-1 bonding restriction imposed 
here was also assumed in Ref . [24[ , which eases compari- 
son of our results to theoretical predictions. 

Two further technical details of the MC algorithm are 
noteworthy: (1) We do not allow any SM pair to both 
break and form a SB (i. e., to break and recombine) dur- 
ing the same MC step. This is a technical violation of 
detailed balance, but satisfies the weaker "balance" con- 
dition sufficient 1811 for accurate MC simulations. SB 
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recombination is a critical aspect of supramolecular poly- 
mer physics [24[ and is further discussed in Section HTH 
(2) we do not allow "bond switching" moves within a 
single MC step. That is, for SMs V, W, A, Y, we do not 
allow moves of the form 



V- w + x 
V -W + X -Y 



V -X + w 
V -Y + X -W 



(2) 



or any other more complicated moves. In addition to be- 
ing difficult to implement in simulations, such processes 
are unlikely to occur instantaneously in real systems, in 
part because of steric constraints. Different but analo- 
gous rules, suitably modified to the use of SMs with two 
binding sites each, were imposed in Refs. [II [77]]. 

In our model, varying the relative rates of sticky bond 
formation and dissociation is accomplished by varying h. 
However, the absolute values of the rates depend on a 
yet unspecified kinetic time scale r^j„ . This time can be 
controlled through the Monte Carlo routine. At each MC 
timestep (i. e., every to), a fraction juc of unbonded SM 
pairs (of those within range r < Rq) and an equal frac- 
tion Jmc of bonded SM pairs are randomly selected to 
be considered respectively for sticky bond formation and 
breaking. We have verified this scheme maintains 'bal- 
ance' [8l| for pairs and triplets of SMs for .01 < Jmc < 1- 

Thus the average time over which each unbonded or 
bonded SM pair is considered once for (respectively) SB 
formation or breaking is tmc = /mc t 0j an( i the param- 
eter tmc effectively controls the "chemical kinetics" of 
the SBs. For a discussion of why we use Jmc < 1 rather 
than varying To, see Section III Dl Small tmc correspond 
to fast chemical kinetics [82| , while large tmc correspond 
to slow chemical kinetics. 

In Section IIII El we perform mechanical tests on vari- 
ous systems. Two types of tests are perfomed; constant 
volume deformation and tensile creep. In the constant 
volume deformation tests, L z is increased at a true strain 
rate e = L z / L z , and L x and L y are adjusted to maintain 
constant volume. In the creep tests, a constant (small) 
stress difference a creep (relative to the equilibrium hydro- 
static pressure in the quiescent state, which is positive for 
repulsive LJ interactions) is applied along the z direction 
using a Nose Hoover barostat [69|. This smaller \a z \ pro- 
duces tensile creep. Both types of tests use r = -2tlj to 
minimize systematic errors. 



B. Comparison to Previous Simulation Protocols 

It is worthwhile to compare the simulation method 
and ambient conditions described above in the context 
of previous AP simulation studies. The use of a hybrid 
MD/MC method is a powerful advantage. Pure Monte 
Carlo (MC) simulations have been performed with lattice 
H3, HE S3, HE HI] and off-lattice [El 53 models. These 
are very effective at studying static properties like per- 
colative gelation and (in the case of solutions) phase sep- 
aration, but have limited ability to capture the complex, 



collective dynamics which occur in bulk polymers, and 
thus lack properties (i) and (ii). For example, MC can 
not, even in principle [39|, capture hydrodynamic effects, 
which are expected to play an important role whenever 
momentum transfer is important (e. g. in relaxation of 
highly stressed systems). 

Pure molecular dynamics (MD) studies fill . l42l . l43l . 
EE EE Hll have been used to study static and dy- 
namic properties. While better able to capture dynamics 
and nonequilibrium phenomena than MC, MD studies 
can not naturally implement controllable sticky bonding 
topology. Also, MD studies cannot easily impose any 
control of chemical kinetics without resorting to costly, 
chemically realistic models. For example, Padding and 
Boek j§6] studied systems intermediate between ours and 
those studied by Huang et. al; a fraction c st < 1 of 
their monomers were allowed to form linear equilibrium 
poiymers, but the FENE-C sticky bonding potential [87| 
used did not allow for variable kinetics. Thus, in practice, 
typical MD studies lack properties (iii) and (iv) . 

The previous works most closely related to the present 
method are Refs. 0, EE Iz3|, w h° also used hybrid 
MC/MD with the same U sb (r) (Eq. [J). Huang et. al. 
[Til Wf\ also used variable kinetics, but studied equil- 
brium linear polymers with c s t = 1 rather than network- 
forming APs with c s t "C 1. Details of the Huang et. al. 
method are discussed extensively in Ref. [lij |. The key 
differences of our method from Ref. [4(| are the imposi- 
tion of binary bonding and the use of variable tmc (they 
used only one tmc = 0.2txj). Another difference was 
that Refs. 0, EElZ2| all used a much stronger thermo- 



stat, giving overdamped (Brownian) dynamics. 

Many previous studies have used nonspecific (e. g. 
strengthened attractive Lennard-Jones or Coulombic) in- 
teractions which allow SMs to form arbit arily many si- 
multaneous sbs [IE EE EE EE EE EE EE M, HI . This 

results in formation of interesting structures such as mi- 
celles and micelle-bridge networks, which occur in real 
AP systems such as associating ionomers (see e. g. Ref. 
[IE]). In contrast, the (experimental) APs we wish to 
model tend to form networks more like classical rubbers. 

Most previous studies [H, HE EE El El, H varied 
temperature T at fixed SM bonding strength. This does 
not isolate the effects of T on associative bonding from 
its other effects such as the dynamical slowdown which 
occurs in normal (non- associative) polymers. To get a 
full picture of AP physics, one should vary both h and T 
independently [l4j]- We follow this approach. 

Other differences from previous simulation studies 
are more associated with the systems employed than 
the methods applied. Many studies have considered 
only telechelic chains 0, El El, El El, EE EE H- 
Telechelics are appealing in their simplicity, but their 
network-forming abilities are naturally limited; for 
binary bonding, at least 3 SMs/chain are required to 
form good networks. Weakly entangled chains (N ~ N e ) 
may be ideal [8| for technological goals such as enhanced 
melt processability at high T and network strength at 
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low T. The majority of previous studies have employed 
extremely short N < N e chains [H E2, El El Efl [Hi , 
but we consider systems with N ~ N e . Finally, the 
majority of previous studies have focused on sm all p 
corresponding to solutions 0, H3, ES El El, HI, HI • 
AP solutions exhibit a wide range of intriguing phe- 
nomena, in particular competition between gelation 
and phase separation [361 |38| |. which, however, we do 
not wish to consider here. In addition, the presence of 
solvent can dramatically weaken the effective strength 
of sticky bonds in real systems [1, [gl, [76j] ; this effect is 
beyond the scope of our model. We therefore focus on 
systems with p corresponding to a dense pure melt with 
no solvent. 

Mappings of the bead-spring model to real, dense 
polymer melts (48[ produce different tlj in the range 
lCP 10 ' 5±1,5 s. Present day computers can achieve runs 
(for the system sizes used here) of up to ~ 10 7 t£j ~ 

10 -3.5±1.5 S) but 

runs this long can not be performed 
over a broad parameter space. In contrast, sticky 
bond lifetimes in experimental systems are typically at 
least 10~fs, and often many orders of magnitude longer 
[Hoi . [60L [63]. Thus any attempt to capture specific 
SM chemistries and at the same time use systems large 
enough to study bulk dynamics and mechanical proper- 
ties would exceed the capabilities of present day super- 
computers [89]. Coarse-grained modelling with the goal 
of studying the dynamics of AP systems by analogy is the 
only currently feasible approach for bulk systems, so we 
make no attempt to mimic specific chemistries. The only 
published simulations of which we are aware that model 
AP networks with specific chemistries [8J, |90| are pure 
Monte Carlo studies that used a very coarse-grained (lat- 
tice) bond-fluctuation model [9l[ and focused on static 
properties. 



C. Static Properties: Validation of Hybrid 
MD/MC Method 

As discussed above, systems contain a total of N st = 
N c hNc s t sticky monomers. Due to the binary bonding 
rules, the maximum number of sticky bonds that can ex- 
ist in the system at any given time is N st /2. If the prob- 
ability that an SM is bound into an SB is p 'active, then 
the total number of SBs in the system is N s tp ac ti ve /2. If 
A represents an unbound SM and A 2 represents a bound 
SM pair, these factors define the concentrations 

[A] = pc st (l -Pactive), 

( 3 ) 

[A 2 ] = pC s tPactive/2, 

where square brackets denote concentrations. If the equi- 
librium value of Pactive is p* , then the equilibrium con- 
stant for SB association is defined (by the law of mass 



action for the reaction A + A <-> A 2 ) as 

\M] P* 
Keq ~ [A]* ~ 2pc at {l-p*) 2 { ) 

for binary bonding. 

Figure [2] shows simulation data in which p* was evalu- 
ated from equilibrated simulations at fixed h and K eq ob- 
tained from Eq. [4j Circles show values of K eq for N = 1 
and N = 50 systems. As expected, K eq ~ exp{h/ksT). 
The data shown are for tmc = 1-Otlj, but we have ver- 
ified that p* is independent of tmc (to within statistical 
errors) for all h tested, over the range tlj < tmc < 
IOOt^j. Because there is an entropy cost ~ ksT to form 
a SB, few SBs form for h < 2u Q . As h/uo ranges from 
2 to 17.5, the equilibrium constant K eq varies over more 
than six orders of magnitude, from 0.96 to 3.9 • 10 6 . This 
is a wider range of h and K eq than considered in previ- 
ous simulation studies. A standard [92] finite size scal- 
ing analysis of the percolation gel transition is given in 
Appendix fA"! For ksT = l.Ouo, percolation occurs at 
h = h perc = 4.25uo, so we consider values of h up to ~ 4 
times above the gelation transition. 

Note that the TMc-hidependence of p* allows systems 
to be equilibrated efficiently using a low tmc = t lj- 
Higher values of h (for polymeric systems) are impossible 
to equilibrate on present-day computers with our current 
method; equilibration is discussed further in Section [Ml 
However, the highest values of K eq considered here are 
comparable to those observed in some experiments on 
multiple- H-bonding SMs [111]. 




FIG. 2: Sticky association in equilibrium; simulation data 
and test of Equation [8] All results are for N ch N = 280000, 
c a t = 0.08 systems with ksT — l.Owo and tmc = I.Otlj. 
Closed circles are simulation values of K eq from Eq. [4] for 
N = 50 polymers and open circles are for TV = 1 dimer- 
forming systems. The straight lines are exponential fits, to 
Eq. U for K^ q s . 

Data from multiple system sizes are also useful in fur- 
ther validating the simulation model. Ben-Nairn and 
Krapivsky have pointed out that systems which rc- 
versibly polymerize undergo a nonthermodynamic gela- 
tion transition [93j when the fragmentation (in our case, 



7 



SB breaking) process is too weak. The average num- 
ber of clusters (aggegrates) at any given time is N agg = 
N c h/N n , where N n is the number-averaged cluster size 
(Appendix N agg = N c h in the absence of sticky 

bonding and N agg — > 1 in the limit of large h, because 
all the chains combine into a single network (as in an 
ideal rubber). Our systems, in the terms of Ref. [93|, are 
"thermodynamic" if and only if: (1) N agg is linearly pro- 
portional to N c h below percolation (i. e. for h < h perc ) 
and (2) the probability distibution of cluster sizes P(M) 
(Appendix IA"]) is independent of tmc- An arbitrary sim- 
ulation method will not necessarily display a 'thermo- 
dynamic' gel transition; failure to do this would be a 
serious flaw according to our goals. We therefore have 
verified that our model satisfies conditions (1) and (2) 
for tlj < tmc < IOOtxj, and therefore properly cap- 
tures reversible gelation. Satisfaction of these conditions 
appears equivalent to the above-verified condition that 
p* is independent of tmc [HI • 

Figure [3] shows data for P(M) at h = 4uq and ksT = 
l.Otto (i. e. just below percolation). The collapse of the 
data shows [93| that cluster formation/dissociation is an 
equilibrium processes and supports our arguments that 
the algorithm satisfies detailed balance for the range of 
tmc considered here. Also, P(M) shows some interesting 
properties which demonstrate that our modelled systems 
form good (rubber-like) networks. The line shows a fit 
to a P oc (Af)~ 5 / 4 power law, which is consistent with 
the fractal dimension Df rac = 4 of aggregrates and the 
expected power law ln(P) ~ — (1 + l/Df rac )ln(M) for 
networks j94|. In contrast, dense telechelic systems have 
an exponential P{M) ~ exp(—Mj < M >) distribution. 
The absence of any large exponential contribution in our 
P(M) at large M indicates that long linear clusters are 
not common. Therefore, though our parent chains only 
contain 4 SMs each, we are confident that that is enough 
to accurately capture AP network physics. 



D. SB Dynamics and Two-State Model 

Figure 0] shows simulation results for the average 
sticky bond lifetime, r s b, in quiescent systems at chem- 
ical equilibrium. Simple thermal activation of SB dis- 
sociation would suggest exponential behavior, t7 oc 
exp(—h/kBT). In fact the results are markedly nonex- 
ponential. Interestingly, SB lifetimes in polymeric sys- 
tems are (apparently) always lower than those in dimer- 
forming systems. This is consistent with differences in 
chain connectivity; SBs embedded in polymers experi- 
ence additional 'pulling' forces due to transmission of the 
random thermal forces (which produce diffusive motion) 
through covalent bonds along their parent chains. Addi- 
tional reductions in r s f, could potentially arise from in- 
creased steric hindrance to bonding for embedded SMs. 
While this "polymeric" effect on r s b should dependent 
sensitively on N, c s t and T, to our knowledge it is not 
included in any theories for AP networks. 




FIG. 3: Cluster size distribution. P(M) is the probability 
that a chain will be part of a disconnected cluster of M chains 
(i. e. the weight fraction of M-clusters). All results are 
for N = 50, uniform c st = 0.08 systems with ksT = ito 
and h = 4uo- Data for different kinetic rates are shown: 
tmc /tlj = 1 (blue stars), 10 (green x), and 100 (red +). 
The upward slope at large MW is due to the statistics of 
small numbers. Results are averaged over 100 statistically 
independent samples. 




h/uo 

FIG. 4: Sticky bond lifetimes. All results are for 280000-bead, 
c 3 t ~ 0.08 systems with tmc = I.Qtlj and ksT — uq. Closed 
circles are simulation data for iV = 50 polymers, open circles 
are data for N = 1 dimer-forming systems, and the straight 
line is an exponential 'fit', shown only as a guide to the eye. 



The simulation data in Figures and 2] can be bet- 
ter understood by mapping the Monte Carlo procedure 
and U s b(r, h) onto a two state Arrhenius model for sticky 
bonding. The model is depicted in Figure O Bonded 
SM pairs are assumed to have an energy —h, unbonded 
SMs have zero energy, and we introduce an /i-dependent 
barrier 5(h). 

The Monte Carlo rules described in Section Hi Al allow 
us to assume that sticky bond formation obeys second or- 
der kinetics and dissociation obeys first order chemical ki- 
netics, as they should as long as pc st « 1/a 3 24, 95] . The 
SB formation/dissociation process can be represented as 
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the chemical reaction 



Refs. QJ, lZl did not intro- 



A + A < — A 2 



(5) 



where kf and kb are the rate constants for SB forma- 
tion and dissociation. Then the equation for chemical 
equilibrium is 



k f [A} 2 = k b [A 2 ] 



(6) 



In the Arrhenius two state model the rate constants are 
given by: 



kf = aexp(—8(h)/kBT), 

k b = 0exp(-(h + 5(h))/k B T), 



(7) 



where a and (3 are constants with dimensions of 
volume x frequency and frequency, respectively. Note 
that the above is a "mean field" model [l2j in that it 
ignores correlations between sticky monomers (i. e. con- 
centration fluctuations). Thus kf and kb (and especially 
a and /3) will in general depend [44| on N, p, c s t, and 
(through second order effects such as the variation of p 
at fixed pressure) T. 

In thermal equilbrium, Eq. Ogives the equilibrium con- 
stant 



K. 



TS 



ki = a 

k b /3 



exp(h/kBT). 



(8) 



Eq. [5] fits simulation results for K eq very well, as shown 
in Fig. [2l In AP networks at even higher values of h, Eq. 
H] should fail due to 'trapped' open SMs [44J that cannot 
find partners, but this effect is negligible for the systems 
considered here. 

We now compare two state model predictions to sim- 
ulation data and map the latter to the former. In 
the two state model, the mean SB lifetime r s & is just 
T sb = k^ 1 . Similarly the probability that an unbonded 
pair in the '2A' state (Fig. O will jump over the barrier 
is just exp(—5(h)/kBT). This is also the success rate 




FIG. 6: Validation of method. All results are for 280000-bead, 
Cst = 0.08 systems with tmc = I.Otlj and ksT — uo. Closed 
circles are simulation data for iV = 50 polymers, open circles 
are data for N = 1 dimer-forming systems, and straight lines 
are exponential fits. Data shown are: (a) S(h)/kBT, and (b) 
r- 1 ex P (5(ft)/fe fl T). 



Smc = kf/a for Monte Carlo SB formation attempts, 
so 5(h) can be directly measured from the simulations: 
S(h)/k B T = -ln(S MC ). 

In Figure panel (a) shows simulation re- 
sults for 6(h) and panel (b) shows simulation re- 
sults for t^ 1 exp(S(h) / ksT) . The latter shows that 
the perfect exponential decay expected from Eq. [71 
T^ 1 exp (5(h) /ksT) = kb = /?exp (— h/ksT), is actually 
observed. Note that this Arrhenius behavior was in no 
way imposed; it emerges naturally, showing the utility of 
the two state model in understanding the behavior of our 
simulations. 

The parameters a and (3 can be extracted from the 
data in Figs. [2] and [6[ For tmc — t l.j, ol — 6.3a 3 /rLj 
and (3 = 24/tlj for dimers, while a = 4.0a 3 /tlj and 
(3 = 41/ tlj for N = 50 chains. The smaller a measured 
for polymer-embedded SMs is consistent with the above- 
hypothesized increased steric constraints. The large rate 
constant [3 = 41/tlj indicates a potential problem with 
the simulations. (3 is an effective "attempt frequency" 
for breaking sticky bonds, which implies that the MC 
timestep To should be small compared to /3 . Larger 
To will in principle produce systematic errors. The data 
shown above are for To = I.Otlj, which is large compared 
to j3~ x . Simply reducing tq is problematic because it 
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sharply reduces the parallel efficiency of the simulations. 

However, we have used values of tq as small as .05tlj, 
and find that all errors produced by using r = 1-0t lj 
are small in quiescent systems at equilibrium; for exam- 
ple, the systematic error in t s {, at h = IOuq is about 1%. 
While the errors in dynamical properties are somewhat 
larger at small h (h -C 10uo), in this paper we focus 
on dynamics for h > lOito and use tq = I.Otlj. In all 
cases, all differences produced by smaller r are small 
compared to the differences between systems contrasted 
in Section IlIIl and comparable to our statistical errors, 
i. e. ~ 1%. For noncquilibrium systems, however, sys- 
tematic errors are larger. Thus all nonequilibrium and 
mechanical-property tests in this paper are performed 
using tq < 0.2t lj . 

In summary, to within our noise, increasing tmc leaves 
the static properties of our model AP networks (Figure 
[2]) unchanged, and changes the sticky bonding dynam- 
ics (Figure [5^-b) only through the prefactor T s b oc T~f c . 
The role of tmc in the dynamics therefore appears in 
the rate constants a and f3, which are also proportional 
to t^ c . Increasing tmc slows down the chemical kinetics 
of the SBs (both formation and dissociation) relative to 
the underlying polymeric time scales, while leaving the 
thermodynamics unchanged. This is why we claim our 
model can separate thermodynamics and kinetics. The 
variation of h and tmc employed here may be thought of 
as corresponding to "scanning" across chemically differ- 
ent sticky monomers. Given the time scale problem men- 
tioned above, this scanning is only qualitative. However, 
we show below that it is very useful in understanding AP 
systems. 



dent on having only a few sticky monomers per chain, al- 
though increasing Nc s t at fixed p will broaden the range 
by lowering h perc . Here we focus on values of h which 
are well below hpcoT, and thus "in the middle" of the 
physical gel regime. 

Another of the key features of physical AP gels is sticky 
bond recombination. The concentration pc s t(l — p*) of 
free SMs is small. Moreover, the motion of free SMs is 
constrained by their (transiently but usually) bonded in- 
trachain neighbors. Thus SM pairs tend to recombinc 
after SB-dissociation events. This leads to to a second 
characteristic timescale for individual sticky bonds; in 
addition to the "bare" lifetime r s 6, there is [24| a larger, 
"effective" SB lifetime t* , which can be thought of as 
the average time for initially bonded SMs to "separate" 
(i. e. no longer recombine) as opposed to merely debond. 
It is of interest because rheological experiments typi cally 
measure r*; T s b is more difficult to access [!, [H, l35j . 
Values for r s b and t* (defined more specifically in Ap- 
pendix [B] and discussed further in Section IIII Bp for a 
wide variety of systems are given in Table |TJ We have 
already shown how T s b is affected by polymer physics 
- indirectly through covalent backbone bonds. Now we 
study the ways in which SB recombination influences and 
is influenced by the interplay of SB thermodynamics, SB 
kinetics, and polymer physics. We perform our study in 
terms of measurements of diffusion, r*, dynamical het- 
erogeneity, nonequilibrium chemical dynamics, and non- 
linear mechanical properties. All results presented below 
are for systems that were first equilibrated for many t*. 
As will be shown, these are best understood by deter- 
mining whether SB recombination is diffusion-limited or 
kinetically limited. 



III. RESULTS 

Previous work has shown (23, Ei| that the most dra- 
matic changes in dynamics, our primary interest, take 
place not at h perc but rather at considerably higher h. 
The rest of this paper considers systems with h 3> h perc 
and p* > 0.95. This is the "physical gel" regime [H[ 
where nearly all chains are (at any moment) part of a 
single aggregate. A "snapshot" of a physical gel looks 
much like a crosslinked rubber, yet chains are delocal- 
ized and the system can flow at long times. One of the 
most interesting properties of physical gels is their tran- 
sition to chemical gels as h increases or T decreases. In 
this "physical-chemical gel transition" (PCGT), chains 
become localized [13, |4|| in a manner analogous to 
the "caging" effe ct p roduced upon cooling fragile glass- 
forming systems [9a H3] • 

For the AT, T, and c st considered here, the PCGT oc- 
curs [98[ at bonding strength hpcGT > 17.5uq. The 
broad range (5uq < h < YJuq) between the percolation 
(Appendix [XJ and localization transitions is consistent 
with the findings of Kumar and Douglas as well as 
Baljon et. al. [46[, who both, however, used constant SB 
strength and varied T . The broad range is not depen- 



A. Diffusion 

The effect of varying different thermodynamic and ki- 
netic parameters on monomer diffusion (mean squared 
displacement < (Sr) 2 (t) >) is shown in Figure [7] Panel 
(a) shows the variation as h is increased at tmc — X.Qtlj 
and ksT = 1 .Ouo- At short times (t -C r s &), results for 
different values of h collapse, showing (as expected) that 
sticky bonding has little effect on diffusion on these time 
scales. At larger times (t > T s b) results show a pro- 
gressive localization and 'caging' effect, similar to that 
described in Refs. [10, E(|, as sticky bond strength in- 
creases. At ft, = 10uo, little localization occurs because 
[2~ij T s b is less than the Rouse time of the chains in the 
absence of sticky bonding {tr — 2.6 • 10 3 T£j). As r s & 
increases with increasing h, the curves develop a "shoul- 
der" which illustrate the temporary caging associated 
with physical gels. This temporary cage becomes per- 
manent as T s b ^ oo (as in a classical crosslinked rubber). 
Data in panel (a) support our earlier statement that this 
occurs for some h > 17uo- 

Panels (b-c) show a pair of interesting effects. First, 
for h = II.25uo and fc^T = I.Ouq, increasing tmc has 
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FIG. 7: Diffusion as function of h, T, and tmc- All systems 
have N c h = 1400. Panel (a): ksT — uo, increasing h. Panel 
(b): h/ksT = 11.25; fcsT = wo and ksT — 0.6uo, increasing 
tmc- Panel (c) is a blowup of (b) showing the crossover. 
Lines from top to bottom for each T in panels (b-c) are for 
tmc/tlj = 1, 10, 100, and (for ksT = l.Otto) oo. The legends 
for panels (b-c) give T in units of uo/ks and tmc in units of 
tlj- 



the same qualitative effect as increasing h at fixed tmc- 
Data for < (5r) 2 (t) > collapse for t less than the small- 
est T s b (i. e. T s b for the lowest tmc)- For longer times, 
the data develops a shoulder which increases in width as 



tmc increases. Data for an equilibrated system with MC 
deactivated (i. e. tmc = 00 ) shows "chemical gel" (ideal 
rubber) behavior; chains are permanently localized. 

Second, data from systems with the same "SB thermo- 
dynamics" (i. e. the ratio h/ksT — 11.25) but different 
ambient conditions (h = 6.75mo and ksT = 0.6uo) shows 
interesting contrasts which illustrate the interplay of SB 
dissociation and underlying polymer physics. For t % T s b, 
data for < (8r) 2 (t) > still collapse, but data from the 
lower-T systems collapse on a lower value. This is not 
at all surprising, as polymeric diffusion is well known to 
slow with decreasing T. However, though h/ksT is the 
same, values of T s b are smaller for the lower-T systems, 
perhaps because h is smaller and SB breaking is favor- 
able at smaller r (see Fig.[I| [93|. Thus the diffusion data 
actually cross over at intermediate time scales (panel c) 
and the lower T systems show greater mobility at fixed 
h/ksT, a most unusual state of affairs. While the case 
presented here is somewhat artificial because in a real 
polymer melt p would decrease with T and lead to further 
diffusive slowdown, we believe the point that varying T 
at fixed SB thermodynamics should change relaxation at 
different timescales differently should be generally valid. 
For example, the freque ncy ( a;) dependence of the dy- 
namical moduli G(uj; T) jl00l | should change with T in 
nontrivial ways. In other words, time-temperature su- 
perposition should be violated. 

It is useful to relate the mean squared displacement 
to the cage size a cage and "escape parameter" f(t) using 
the definition 



< (tf0 2 (*) 



cage 



fit). 



(9) 



In the T s b —>■ oo limit, o? cage y (p ent) is the volume 
explored by sticky monomers jl0l| . The "chemical gel" 
time t c hem is the time at which f(t) approaches unity 
in this limit. Data for the tmc = oo system in Figure 
\3b) (with h = 11.25u , N = 50, and c st = .08) shows 
that a 2 cage ~ 23a 2 and t chem — \Q A b T LJ . For finite T sb , 
one can define t cage as a "caging" time describing the 
(de)localization of SMs [4(| EH; f(t) then has the general 
form /(0) - 0, /(<) ~ 1 for t ~ t cage , and /(*) > 1 for 

t ^ tcage 

The monomeric diffusion constant D, as measured 
by limt^oo < (5r) 2 (t) >~ 6Dt, should vary inversely 
with some "long" characteristic time Ti ong of the sys- 
tem, roughly defined as the time for chains to diffuse by 
their end-end distance. Candidates for Ti ong include T s b 
and t*. Ref. [2~i| predicts Ti ong oc T s b for weakly bind- 
ing physical gels and Ti ong oc r* in the strong-binding 
(near-chemical) limit. Figure [8] shows results for DT s b, 
DT star , and 8 • 10 4 r L ,/D from Table [fl for h = 11.25w , 
ksT = l.Ouo systems, over a wide range of tmc- D de- 
creases with increasing tmc slower than both (t*) _1 and 
T~ b , but it tracks the former more closely than the latter. 
Thus results for this value of h are apparently intermedi- 
ate between the "weak" and "strong" physical gel limits 
described in Ref. 12411. 
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FIG. 8: Scaling of diffusion with various candidate "long" 
relaxation times. Data shows Dri ong for r; ong = T s f, (circles), 
Ti orlg = r* (triangles), and Ti ong = 8 ■ 10 4 t L j (squares). The 
last value is chosen so the "bare" diffusion constant D can 
be shown on the same plot. All results are for systems with 
N ch = 1400, h = 11.25u and k B T = l.Ouo- Dtmc is not 
shown because T sb oc tmc- 



B. Crossover between Diffusion-Limited and 
Kinetially Limited SB recombination 

Table Q] shows r* and r* /r s & for all investigated sys- 
tems. As expected, increasing /i at fixed tmc increases 
T * / T sb, because for fixed kinetics recombination is more 
likely for thermodynamically stronger SBs. There are 
several possible regimes of possible relations between r* 
and r s b that can be related to diffusion (specifically, 
< (8r) 2 (t) >) on intermediate timescales. Systems with 
T s b 3> fnin{tcagejt c hem) wlu exhibit kinetically limited 
sticky bond recombination (KL); in this regime r*/r s b is 
predicted to be constant [24(. However, if r s b -C r cage , 
recombination will be dominated by "correlated" recom- 
binations of SM pairs that have recently dissociated and 
have not had time to fully diffuse away from one an- 
other. This is diffusion limited sticky bond recombina- 
tion (DL). To our knowledge, the DL regime and espe- 
cially the crossover between DL and KL have not been 
previously studied for AP networks. 

Figure [5] shows the variation of t*/t s (, with chemical 
kinetics for h = 11.25uo systems at ksT — I.Ouq. Over 
a range of two orders of magnitude in t mCi the data are 
well fit by the equation 
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Tsb 



C 



K 



(jMCr 



(10) 



where C is the probability of recombination in the KL 
limit and K is the contribution from diffusion-limited re- 
combination. The exact form of Eq.[l0]is not of great con- 
sequence; what matters is the broad crossover between 
regimes and the large change of r*/r s b as a function of 
kinetics. Nevertheless, since T sb = fcr 1 , Eq. \TU\ can be 
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FIG. 9: Crossover from diffusion limited to kinetically limited 
sticky bond recombination for h — 11.25. All systems have 
70000 beads, with c st = .08. The upper data set is for N = 50 
polymers and the lower is for N — 1 dimer-forming systems. 
The fastest-kinetic systems (tmc ~ 0.5tlj) use to — 0.5tlj, 
and all others use to — tlj- The solid line shows a fit to Eq. 
HO] with C = 1.15, K = 4.98rfj and x = 0.74. The dashed 
line shows a fit with C set to 1, K = 2.97rf j and x = .94. 



interestingly rewritten 



C 



Kt~ x 



(11) 



The significance of Eq. [TT] is its prediction of a nonlinear 
dependence of kbT* on the rate constant fcf, for dissocia- 
tion; t^ c cx kf (see Section Hi Dp . 

C and K are of course not universal constants, but 
will depend on h, p, c stl T, and N. In practice, one 
would expect Kt% < C when T sb > min(t cage ,t chem ). 
More physically, the condition Kr^ b <C C defines the 
KL regime, where sticky bond reactions become "mean- 
field" in the sense of Cates [l2[ • The data in Figure [9] 
show that one can (at least in our model systems) move 
from the KL to the DL regimes simply by speeding up 
the chemical kinetics, if one is in the regime where T s b is 
comparable to the underlying polymer relaxation times 
such as t c f lem . 

For our systems, values of C are close to values of 
t*/t s i, in our "kinetically slow" (tmc — IOOtxj) systems, 
cf. Table HI C increases with h (qualitatively) as pre- 
dicted by Rubinstein and Semenov [24[. However, while 
our kinetically slow systems all have C < 2, the "strong 
physical gel" theory in Ref. [24| assumes C ^> 1, so we 
defer a detailed comparison to that theory to later work. 
Here we merely make the positive observation that the 
basic prediction [24j of effective SB lifetime renormaliza- 
tion (r s h — * t* ) works well at these relatively small C and 
(somewhat surprisingly) over the entire studied KL— >DL 
crossover regime. The renormalization r s b — * r* accu- 
rately captures e ffectiv e SB dissociation over a very broad 
parameter space [l02j |. 

On the other hand, an observation apparent from Fig. 
[5]is that SB recombination is in general only partly 'poly- 
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meric' in nature. In the limit of fast kinetics, values 
for t*/t in N = 1 systems are nearly as high as for 
N = 50 systems. However, the lack of a network pre- 
vents any caging effects, and so r* /r s 6 decreases much 
faster with increasing tmc, approaching unity (regard- 
less of the value of h) at tmc = IOOtlj. This illus- 
trates that the crossover from DL to KL is analogous to 
a crossover between 'dimeric' and 'polymeric' recombina- 
tion; in other words, chain connectivity (i. e. covalent 
bonding) becomes increasingly important as kinetics are 
slowed. While the 'dimeric' contribution to SB recom- 
bination cannot be simply "subtracted out" due to the 
different 6(h), these 'dimeric- vs-polymeric' effects on SB 
recom bination have been neglected by previous theories 

Ho3. 



C. Recombination and Dynamical Heterogeneity 




FIG. 10: SB recombination and dynamical heterogeneity. All 
systems have N c h = 1400, h — 11.25uo and ksT = l.Ouo- 
Lines from top to bottom are for tmc/tlj = 1, 10, and 100. 
Panel (a): P r ecomb(t). The black dotted line is exp(—t/r*) — 
exp(—t/T s b) for tmc — Wtlj. Panel (b): A(t). 

Figure [TUh shows simulation results for the SB recom- 
bination probability P re comb(t) (defined in Appendix |B|) 



for h = 11.25uo systems with different chemical kinetics. 
Although there is a small peak in P reC omb (not displayed 
and low compared to the peaks shown in Fig. HDh ) at 
very small times t ~ tmc, the i~ 5 / 4 behavior predicted 
[ 1 31 ] for the extreme diffusion-limited case is not found. 
This indicates none of our systems have "too-fast" kinet- 
ics [l3|. For t ^> To, our results have the interesting form 
Precomb(t) — exp(— t / t*) — exp(— t/T s b)', a comparison to 
actual data for tmc = IOtxj is shown. 

This form of P re comb(t) has a maximum at the "der- 
ealization" time 



T s b- 



yiogj/ 

1 ' 



y 



(12) 



where y = t* /r s b. Bonded SM pairs trend towards mov- 
ing away from each other after Tdeioc- Furthermore, 



7~delc 



logy 



y 



r 



(13) 



which becomes small for y ^> 1. This suggests Tdeioc 
(in addition to t*) might be a key relaxation time in 
systems with large y. However, this is speculative and 
needs further verification. 

A useful measure of relaxation in complex fluids is the 
"non- Gaussian" parameter 



A(t) = Z< ^> 2 - 1, 
w 5 < 8r 2 (t) > 2 



(14) 



which is zero for normal diffusion and positive for sys- 
tems where some particles move anomalously fast [42j . 
particularly for "hopping" type motion. A has been 
shown to be relevant to the structural relaxation of super- 
cooled liquids and dynamical heterogenity [96] . The time 
tdeioc at which A is maximized and the maximum value 
A raol = A(tdeioc) both increase with decreasing T in var- 
ious systems, including associating polymers [40. l42l. |4(| . 
as localization increases, tdeioc may be regarded as a 
crossover time after which the system begins to show liq- 
uidlike behavior. 

Figure [T0b shows the effect of kinetics on A(t). The 
effect of slowing kinetics at fixed h/ksT is similar to the 
effect of increasing h/ksT observed in previous studies 
[40l HfiJ . It is interesting that increasing tmc increases 
dynamical heterogeneity. The probable reason is that 
increasing tmc? even though it leaves p* unaffected, de- 
creases the likelihood of multiple closed SBs on the same 
chain breaking within a short time period. This is consis- 
tent with the idea [2l[ that coherent breaking of nearby 
SBs along a chain eases large-scale motion. The increas- 
ing dynamical heterogeneity with increasing t cage is con- 
sistent with other results showing A max increases as lo- 
calization "transitions" are approached, e. g. stretched- 
exponential relaxation of finite clusters [4lT ]. 

The data in Figs. [51TTU1 also clearly show that tdeioc — 
Tdeioc in systems where recombination is likely (i. e. when 
r*/r s fc is large compared to 1), and thus that delocaliza- 
tion is closely related to individual sticky bonds find- 
ing new partners in a "hopping" type motion. However, 
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these delocalization times are large compared to t cage . 
This is not surprising, as full delocalization should occur 
only when chains have lost all memory of their initial SB 
topology; this "memory" time is inherently polymeric in 
that it must increase with increasing Nc s t, similarly to 
a Rouse or reptation time [4|. Interestingly, the peaks 
of A are broader than those of P re comb- This also likely 
arises either from cluster effects [24| EH or other underly- 
ing many-SM phenomena that ultimately arise from the 
'polymer physics', i. e. the covalent connectivity of the 
parent chains. 



D. Nonequilibrium Chemical Dynamics 

An important feature of our model is its ability to ac- 
curately capture the dynamics of systems in which the 
sticky bonds are not in thermal equilibrium. The evolu- 
tion of SB concentration is, following Equation [5[ given 
by 

[A 2 ] = k f [A] 2 - h[A 2 }, (15) 
which after plugging into Eq. [3] and simplifying becomes 
P active = 2kfpc st (l P active 

r - hp active • (16) 

Equation [TTj] has an analytic solution. For the special 
initial condition p a ctive{t) = at t = 0, the solution is 



Pactive(t) — d 



(d 2 - 1) tanh (2zVd 2 ~^t) + d^/d 2 - 1 
dtanh (2z^d 2 - It) + Vd 2 - 1 



(17) 

where z = pc s tkf and d = 1 + fcf,/4z. 

Figure [TT1 compares this analytic prediction to simula- 
tion results for p a ctive(t) upon activation of sticky bond- 
ing for two systems with the same value of h/ksT but 
different values of h and ksT. Values of z and d in Eq.fTTl 
are taken from fit values of a, (3 and the measured value 
of 8(h) as reported in Section [TTJ note that these vary 
somewhat with T, giving different p* at the same h/ksT. 
Data agree excellently with predictions at short and long 
times. The merely qualitative agreement at intermedi- 
ate times is no cause for concern, but is an interesting 
'feature', because Eq. [TBI ignores all physics arising from 
the important fact that the sticky monomers are embed- 
ded, at a concentration c s t, in chains of length N, in a 
dense polymer melt. The slower convergence of simu- 
lation results for p ac tive(t) relative to the prediction of 
Eq. [T7] is consistent with such polymeric effects; better 
agreement is observed for dimer systems. As expected, 
the polymeric slowdown is greater at lower T. 

The results above demonstrate the ability of our 
method to capture the effect of polymer physics on 
nonequilibrium "chemical dynamics". Thus, as in Ref. 
[T3 |. it can be used to perform "T-jump" simulations. 
These may be useful in analyzing phenomena observed 
in recent real T-jump experiments; nonequilbrium sticky 
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FIG. 11: Nonequilibrium capability of method: /i-jump. Solid 
(dashed) lines are the predictions of Eq.[T7]and upper (lower) 
circles show simulation data for h = 10uo, ksT = l.Owo (h = 
6mo, ksT — 0.6tio) for an N c h = 5600 system after turning 
on sticky bonds. Simulations used to = .05tlj- 



bond behavior is also expected to play a role in self heal- 
ing AP systems @. Note, for example, that the timescale 
over which p ac uve changes in Figure [TTJ is smaller than 
the equilibrium T s b (1.5 • 10 3 tl,/ for h — 10uq). Similarly, 
the timescale of self healing at a fractured surface (where 
Pactive is out of equilibrium) was found to be smaller than 
the time scale for near-equilibrium creep relaxation Q . 



E. Nonlinear and Nonequilibrium Mechanical 
Properties 

Figure [T!2a shows results for creep tests of two sys- 
tems with different thermodyamics and kinetics but the 
same T s b (r s b ~ 1.5 • 10 4 Tr,j). Both tests were per- 
formed at kftT — l.Ouo- The applied stress difference 
|o"z — {?x + cr i/)/2| = .01uo/a 3 is small. At times t -C T s b, 
the extension ratio X z = L z /L a z is the same for both 
systems. For t 3> r s 6, X z is nearly linear in ln(t), im- 
plying that the flow is nearly-linear creep. The system 
with stronger bonds and greater SB recombination shows 
greater resistance to flow, i. e. a smaller creep compliance. 

It is interesting to relate the creep response to the qui- 
escent dynamics. Figure [T^b shows (quiescent) diffusion 
in the same systems. The creep response and diffusion 
are remarkably similar; the onset of more rapid creep in 
the h = lOuo, t mc — IOtlj system under stress corre- 
sponds directly to the onset of (relative) delocalization 
in the quiescent state. This is consistent with a recent 
experiment showing connections between creep behavior 
and linear rhcology in reversible supramolecular networks 

1- 

Next we consider constant volume tension simulations 
at h = 11.25wo, fcsT = uq, and various tmc- These 
simulations can be considered to be an extension of Ref. 



104| , which allowed breaking and formation of interchain 



bonds only at a few (discrete) strains; here SBs break 
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FIG. 12: Linear creep and quiescent diffusion for two systems 
with the same r a b but different SB recombination. Results 
are for systems with N c h = 5600, N = 50, and ksT = l.Owo- 
Panel (a) Stretch \ z — L z jL% under a creep stress Act = 
.01wo/a 3 . The upper curve is for h = 10uo, tmc = Wtlj 
and the lower curve is for h = 12.5uo, tmc ~ I.Otlj. Panel 
(b): mean squared displacement for the same systems in the 
quiescent state. 



standard nonlinear rubber elasticity 105, 106]). 
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FIG. 13: Constant volume tension simulations for systems 
with different SB kinetics. The stress difference Su — a z — 
{a x + cr y )/2 is plotted against <;(A) = A 2 — 1/A. Systems 
have N ch = 5600, N = 50, h = 11.25u , k B T = 1.0u , and 
the strain rate is e = A/A = 10~ ' /tlj. These runs use 
ro = .2tlj for greater accuracy. Data from top to bottom 
correspond to tmc = oo (no SB breaking/forming allowed 
during deformation), tmc = IOOtlj and tmc = 10 Tlj- 
The solid lines are predictions of Eq. 1181 with the value of G e 
taken from a fit to the tmc = oo data and values of r s b taken 
from Table |U 

The simplest result, assuming that e _1 3> Tf ree , 
T s b ^ Tf ree (here Tf ree is the lifetime of unbonded SMs), 
and that SB breaking and formation rates do not vary 
with stress /strain, so that stress memory is lost like 
exp(-t/r s b), is [HI 



er z (A) = G e exp(-ln(\) /er sb )x 
1-A 2 1/A-l 



5(A) 



1 - 2er, 



si, 



1 + £T sb 



(18) 



and reform continuously. Here we present results for 
e = 10~ 5 ' 5 /tlj. Simulations at other e were considered; 
larger values e > t~^ ge make the non-SB-related viscous 
stress contribution unacceptably large, while smaller val- 
ues lead to more sticky bond breaking/formation during 
deformation than is desirable at the values of h and tmc 
considered. 

Figure [13] shows the stress difference | a z — (a x + cr y ) /2 1 . 
With MC deactivated during deformation (fuc — or 
equivalently tmc = °°)j the stress ta kes a form close to 
that predicted by entropic elasticity jl05| : a = G e g(X), 
where g(X) = A 2 — 1/A and A = L Z /LP Z as above. G e 
is predicted to be NcstPinterksT '/2, where p mter < P* 
is the interchain portion of active SBs; the actual value 
from the fit, G e = .028uo/a 3 , is close to the predicted 
value .03lMo/a 3 , indicating viscous stresses are low at 
this strain rate. The fit is performed for g < 5; the 
nonlinear behavior observed at higher g arises from finite 
extensibility of chain segments between crosslinks (as in 



where the first term in brackets is classical rubber elas- 
ticity and the second two terms reflect new SBs created 
during deformation. The Zn(A) comes from the constant 
true strain rate e = dln(X)/dt. 

In Figure fT3l stress-strain results from simulations are 
compared to predictions from Eq. [18] using values for r s b 
from Table [J and the value of G e from the tmc = 00 
system are shown. We confine the comparison to the 
linear regime (g < 5) to avoid confusion. Sticky bond re- 
combination might be expected to slow relaxation. How- 
ever, stress relaxation is actually faster than predicted 
by Eq. [181 We have verified that p ac tive does not de- 
crease during deformation. It appears that instead, T s b 
is reduced by stress. A detailed examination of this ef- 
fect and comparison of nonlinear mechanical properties 
to theories, e. g. Ref. [28| and transient network models, 
e. g. Refs. [Hjlil, HH, is deferred to later work, but the 
data presented above suggests traditional theories will 
break down in the nonlinear regime. 
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TABLE I: Variation of Tab, t* and D with TV, T, h and tmc- Times are in units of tlj. Statistical errors are roughly ±2% or 
less. All systems have N c h = 1400 and data are averaged over multiple statistically independent states. * denotes to = .5tlj 

results. indicates calculation is prohibitive or we have insufficient data. Results for D in TV = 1 systems are not presented 

because they are negligibly affected by sticky bonding. 



System 


TV 


k B T/u 


h/k B T 


tmc 


Tsb 


r* 


T*/T sb 


W 5 T LJ D/a 2 


A 


50 


1.0 


10 


1 


1.55 ■ 10 3 


7.05 ■ 10 3 


4.5 


18.3 


B 


50 


1.0 


10 


10 


1.54 • 10 4 


2.63 • 10 4 


1.7 


8.44 


C 


50 


1.0 


10 


100 


1.53 • 10 5 


1.78 • 10 5 


1.2 


2.17 


E* 


50 


1.0 


11.25 


0.5 


2.40 ■ 10 3 


2.28 ■ 10 4 


9.5 


9.31 


DE* 


1 


1.0 


11.25 


0.5* 


3.35 ■ 10 3 


2.23 • 10 4 


6.7 


N/A 


F 


50 


1.0 


11.25 


1 


4.82 ■ 10 3 


2.95 ■ 10 4 


6.1 


6.89 


DF 


1 


1.0 


11.25 


1 


6.69 • 10 3 


2.67 • 10 4 


4.0 


N/A 


G 


50 


1.0 


11.25 


2 


9.78 ■ 10 3 


4.05 ■ 10 4 


4.1 


6.07 


DG* 


1 


1.0 


11.25 


2 


1.34 ■ 10 4 


3.46 ■ 10 4 


2.6 


N/A 


H 


50 


1.0 


11.25 


5 


2.40 ■ 10 4 


6.42 • 10 4 


2.7 


4.81 


DH* 


1 


1.0 


11.25 


5 


3.36 ■ 10 4 


5.50 • 10 4 


1.6 


N/A 


I 


50 


1.0 


11.25 


10 


4.73 ■ 10 4 


9.70 ■ 10 4 


2.1 


3.11 


DI* 


1 


1.0 


11.25 


10 


6.42 ■ 10 4 


8.52 ■ 10 4 


1.3 


N/A 


J 


50 


1.0 


11.25 


10 15 


1.55 ■ 10 5 


2.32 ■ 10 5 


1.5 


1.95 


DJ* 


1 


1.0 


11.25 


10 1 ' 5 


2.05 ■ 10 5 


2.30 • 10 5 


1.1 


N/A 


K 


50 


1.0 


11.25 


100 


4.82 ■ 10 5 


6.28 ■ 10 5 


1.3 


0.94 


DK* 


1 


1.0 


11.25 


100 


6.70 • 10 s 


6.85 ■ 10 s 


~ 1.02 


N/A 


L 


50 


0.6 


11.25 


1 


1.88 ■ 10 3 


8.02 • 10 3 


4.3 




M 


50 


0.6 


11.25 


10 


1.91 • 10 4 


3.15 • 10 4 


1.6 




N 


50 


0.6 


11.25 


100 


1.93 ■ 10 5 


2.22 ■ 10 5 


1.2 







50 


1.0 


12.5 


1 


1.53 ■ 10 4 


1.23 ■ 10 s 


8.0 


2.79 


P 


50 


1.0 


12.5 


10 


1.55 ■ 10 5 


3.58 ■ 10 5 


2.3 




Q 


50 


1.0 


12.5 


100 


1.53 • 10 6 


2.07 • 10 6 


1.3 




R 


50 


1.0 


13.75 


1 


4.90 ■ 10 4 


5.42 ■ 10 s 


11 




S 


50 


1.0 


13.75 


100 


4.62 ■ 10 6 


7.34 ■ 10 6 


1.6 




T 


50 


1.0 


15 


1 


1.56 ■ 10 5 


2.49 • 10 6 


16 




U 


50 


1.0 


15 


100 


1.48 • 10 7 


2.53- 10 7 


1.7 




V 


50 


1.0 


16.25 


1 


5.0 ■ 10 5 


1.1 • 10 7 


22 





IV. DISCUSSION AND CONCLUSIONS 

We have performed an initial set of simulations us- 
ing a new coarse-grained model for associating polymers. 
The MD /MC hybrid algorithm and variable chemical ki- 
netics allow for greater realism and flexibility than in 
previous simulations of AP networks. Further, the 1-1 
sticky monomer binding topology imposed here reflects 
current experimental trends. The model was extensively 
validated and is able to accurately model equilibrium 
dynamical properties, nonlinear mechanical properties, 
and far-from-equilibrium systems. We studied the model 
over a very broad parameter space. While have empha- 
sized that we study APs by analogy because simulations 
of chemically realistic AP networks are not yet compu- 
tationally feasible |89l|. our results should nevertheless 
aid in "rational" [5l| design of AP systems, especially in 
"transition" regimes like those discussed in this paper. 

The key results presented here focused on separa- 
tion, comparison and contrast of thermodynamic and 
chemical-kinetic effects on SB recombination, the mo- 
tion of individual chains, and bulk mechanical proper- 
ties. As expected, instantaneous network structure was 



independent of kinetics at fixed thermodynamic condi- 
tions (i. e. sticky bond strength /i/fcgT), and relaxation 
times increases with increasing h/ksT. Similarly, at 
fixed h and IcbT, relaxation slows as the chemical kinet- 
ics are slowed. This was illustrated by measurements of 
monomer diffusion. In the physical gel regime, monomers 
experience a temporary "caging" similar to that found in 
glasses. This caging effect strengthens as SB strength is 
increased 0, Sal but also as kinetics are slowed at fixed 
SB strength. Analyses showed that chains become in- 
creasingly localized in a manner similar to that associated 
with the increase in dynamical heterogeneity in non-AP 
melts approaching the glass transition. Of course, the 
analogy should not be taken too far; in AP networks 
the caging is produced only by sticky monomers while 
in systems approaching T g it is produced by hard core 
repulsions of all monomers. 

We find, as expected, that the chemical kinetics con- 
trolling r s t are "mean-field" [ijj and mappable to a two- 
state Arrhcnius model. However, as kinetic rates are in- 
creased and SB recombination becomes non-kinetically- 
limited, the relation between the SB lifetime r s b and 
other relaxation times, such as the effective SB lifetime 
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r* , becomes decidedly nontrivial. This was explicitly re- 
lated to the crossover to diffusion limited SB recombina- 
tion. A new quantitative relation between r* and r s b was 
found. Such relations should be of interest because rhe- 
ological experiments can typically only access t* , which 
is assumed to control stress relaxation (e. g., because 
scission followed by quick recombination does not relax 
stress) dill HI]. 

While the results for T- variation and mechanical prop- 
erties presented here were limited and somewhat prelim- 
inary, we showed examples which illustrate important ef- 
fects. Analysis of diffusion on intermediate time scales il- 
lustrated the point that sticky bond and underlying poly- 
meric timescales will in general vary differently with T, 
affecting the 'interplay' in nontrivial ways. In two sys- 
tems with the same t s {,, systems with different propensi- 
ties for SB recombination showed the same creep flow at 
short times, but those with greater recombination showed 
a smaller long-time ("DC") creep compliance. These dif- 
ferences were directly related to the faster derealiza- 
tion of chains in the quiescent state for systems with 
less recombination. Constant-volume deformation stud- 
ies showed that, as expected, r s (, is reduced by stress. 
Extensive studies of the variation with T (in systems in- 
cluding attractive nonbond interactions for greater real- 
ism) and more detailed analyses of nonlinear mechanical 
properties are underway. 

Nearly all published analytic theories for AP networks 
assume a single controlling relaxation time, either r s & 
or t* , controls the ultimate relaxation properties (i. e. 
other relaxation times scale with the controlling time). 
We showed though various measurements that there is 
a broad parameter space (both in terms of SB strength 
and kinetics) within the physical gel regime where the 
"scaling" assumption fails. This parameter space cor- 
responds to the conditions (1) r*/r s 6 is larger but not 
"much larger" than unity, and/or (2) SB recombination 
is not kinetically limited. Deviations from this "scal- 
ing" behavior due to multiple controlling relaxation times 
have been observed [5l|, [58| , but had not yet been well 
understood. These deviations had been previously as- 
sumed to arise from chemical disorder, and this is no 
doubt partially correct, but as discussed in this paper, 
they also arise from the 'interplay' between SB thermo- 
dynamics, kinetics, and polymer physics. If either (1) or 
(2) hold, both traditional [12, H|| and more sophisticated 
pill 12-il |26| theories should fail to predict the mechanical 
properties of AP networks. This is not meant as a crit- 
icism of the theories, merely an observation that there 
is a broad parameter space where one or more of their 
assumptions fail. 

The DL and KL limits have been discussed by 
O'Shaughnessy and Yu [13| : they respectively correspond 
to dominance of the i-T-term and C-term in our Eq. (TTUJ). 
Conditions under which systems may lie outside the KL 
limit and/or evidence for systems which lie outside it 
are also discussed, to some extent, in the context of AP 
networks in Refs. [H, 0, Coupling between SB 



and polymeric relaxation has also been treated approx- 
imately by Cates [Hj] and Leibler et. al. [2l[, respec- 
tively for linear EP systems and AP networks where 
recombination is improbable. Among published ana- 
lytic theories for AP networks, Refs. [Hi [H| qualita- 
tively treat non-kinetically-limited systems and Ref. [13] 
treats SB recombination. Our results are consistent with 
the argument of Ref. [55| that reaction rates (here de- 
fined as non-recombinative SB exchange) reach the mean- 
field/KL regime only when reaction is slow compared 
to the longest "underlying" p olym eric relaxation time 
[in our case min (t cag e 7 t c hem)\ [107] . Interestingly, Refs. 
[52l l54j suggest that MF kinetics would apply to r* as 
well as T s b in dimensions d > 4 because < |<5F| 2 >^ t y 
necesarily has y < 4 (i. e. no diffusion-limited regime is 
possible) . This suggests that diffusion- limited SB recom- 
bination, which increases t*/t s 6, will increase in impor- 
tance in AP systems with effectively reduced dimension- 
ality (e. g. very thin films or "pores"). Combining the 
approaches of Refs. [U [24], [55[ may be useful for devel- 
oping optimal analytic theories of these systems, at least 
for T well above T g . However, we are not aware of any 
quantitative discussion of the crossover betw een the DL 
and KL regimes such as presented here [l08| |. 

The rheologically simple (i. e. all key relaxation times 
scale with r s b (50ll5l|) behavior observed in the majority 
of experiments on AP networks indicates they exhibit KL 
behavior. Note that these experiments have shown KL 
behavior even though their values of p* are comparable 
to values for systems which in our model exhibit DL be- 
havior at low tmc ■ This likely arises from the slow kinet- 
ics caused by the bulkiness and directional interactions 
of real sticky monomers. Creating (real) strong-binding 
SMs with even faster kinetics seems to be difficult. How- 
ever, one can move out of the KL regime (at fixed kinetic 
rates) simply by slowing the polymeric relaxation times, 
e. g. by going to higher concentrations and/or entangled 
chains. Experiments in this regime are underway [58|, 
and seem to show a breakdown of the simple scaling; 
for example, they show an unusually high power law de- 
pendence of viscosity on concentration, which appears to 
arise because r,;, is of order the time scale for reptation. 
Ref. [5l| also shows an apparent (if weak) breakdown in 
scaling at the highest frequencies considered. In the con- 
text of these observations, we note that the parameter 
space where (1) and or (2) hold may be of greatest inter- 
est for designing materials with novel mechanical proper- 
ties. Our model seems well suited to aid in understanding 
the complicated behavior of AP networks in this regime. 

Here we have left the regime of physically entangled 
APs, which is the re gim e treated by some key analytic 
AP theories [U, [H, [23], untouched. Also, the "inter- 
play" described in this paper should depend on the de- 
tails of sticky monomer arrangement along chains, not 
just N and c s t- Studies of systems with a wide range of 
N and c s t, as well as inhomogeneous (chemically disor- 
dered) systems, are underway. 

In real AP networks the sticky and regular monomers 
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have different sizes and chemistries. Thus an obvious ex- 
tension of our model would be to increase the differences 
between the sticky and normal monomers. For example, 
changing seco ndary interactions may induce microphase 
separation 0, 109]. Another extension would be use of a 
more realistic sticky bonding pote ntia l such as those used 
to model H-bonds (see e. g. Ref. [llOj and refs. therein), 
but here we have focused on chemistry-independent prop- 
erties. A coarse-grained way to to capture this would be 
to keep the binary bonding rules, but include more than 
one sticky site per SM; this would increase the direction- 
ality of bonding, which is a key to the performance of the 
UPy syste ms I 6CH] . Finally, nanocomposites of associating 
polymers 11 1 ill . wher e the presence of nanoparticles may 
or may 

not [ll| [HI 

affect single-chain structure but will 
certainly affect AP network structure, should have even 
richer physics than regular polymer nanocomposites. 
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APPENDIX A: PERCOLATION GEL 
TRANSITION 

Let Nf = N c / l P(N c i l , i) be the average number of dis- 
connected clusters of i chains in a system of N c h chains 
with periodic boundary conditions (at any given time). 
P{N c h,i) is a cluster size probability distribution with 
iP{N c h,i) = 1. The number averaged cluster size 
is then N n — N c h J2i=i iP(N c h, i) and the weight av- 
eraged cluster size is N w — N 7 ^ 1 N ch Yli=i i 2 P(N c h, i)- 
For an infinite system, the percolation gel transition oc- 
curs (by definition) when p* exceeds p perc ; N w diverges 
at p* = Pp erc [49J. However, computer simulations are 
limited to finite N c h, and the value of N w cannot ex- 
ceed N c h- Thus the A^/j-dependent geometric percolation 
p* = Pspan (at which one aggegrate spans the system) 
approaches p per c from below as N c h — > oo) [Hj]. Fortu- 
nately, Pp erc and h perc (the value of h at which p* = p per c 
in an infinite system) can be estimated for finite N c h us- 
ing a standard finite size analysis 69]. 

We perform such an analysis, following Ref. [§2]. Fig- 
ure H3] shows this analysis for N = 50, c st = .08 sys- 
tems at ksT — l.Ouo- The figure plots the rescaled 



variables (N w /N ch )-> N^ 1 ^ vs. (( Pperc - p*)/p*) N^ 3u 
[92 ]. pP erc ~ .40 is close to the predicted value 
l/(Nc st p*f ln ter - 1) [II, where f mter is the fraction of 
SBs which are interchain rather than intrachain. The ex- 
ponents used to collapse the data in the figure are 7 ~ 1.7 
and v ~ 1.2; considering a narrower range of h gives val- 
ues consistent with predictions from the theory of critical 
phenomena (7 ~ 1.8, v ~ 0.9) [l!4| . These ex pon ents 
been extensively discussed in the literature [38l 164 l92j 
and need not be discussed further here. In the figure, h 
increases going from right to left, and percolation occurs 



at h 



P erc 



4.25u - 




((Pperc-P*)/Pperc)N C h V 

FIG. 14: Finite-size analysis of the percolation gel transition. 
Data are at k B T = l.Owo for N ch = 700 (circles), N ch = 1400 
(squares), N c h = 2800 (triangles), and N c h — 5600 (dia- 
monds) are shown. p pe rc ~ .40. 



APPENDIX B: NUMERICAL ANALYSIS OF SB 
RECOMBINATION 

In this paper, the sticky bond self-correlation function 
P altto (A£) is the probability that a bond between two 
given SMs exists both at times t and t + At, while the 
SB "transition function" P trans (Ai) is the probability 
that the bond exists continuously between times t and 
t + At. The SB "recombination function" P r ecomb(At) = 
P au t (At) — Ptrans(At) is the probability that a pair of 
SMs will be bonded at two times separated by At but 
that the bond between them has broken at least once dur- 
ing that interval. We find Ptrans exhibits nearly single- 
exponential decay, Ptrans = exp(— t/r^b) for all systems, 
while for h > 10uo, P a uto also shows exponential decay 
for (at least) the first decade. Values of r* presented 
here are measured from fits to Ptrans = exp(— t/r*). All 
quantities are averaged over all SM pairs and all t. 
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